inverse Subroutine

private subroutine inverse(gamma, nest, inv_gamma)

Inverse matrix Method: Based on Doolittle LU factorization for Ax=b

input: gamma(nest,nest) - array of covariance coefficients for matrix gamma nest - dimension outpu: inv_gamma(nest,nest) - inverse matrix of A

IMPORTANT: the original matrix gamma(nest,nest) is destroyed during the calculation. Define A for avoiding this.

Arguments

Type IntentOptional Attributes Name
double precision, intent(in), DIMENSION(nest,nest) :: gamma
integer, intent(in) :: nest
double precision, intent(out), DIMENSION(nest,nest) :: inv_gamma

Variables

Type Visibility Attributes Name Initial
double precision, public, DIMENSION(nest,nest) :: A
double precision, public, DIMENSION(nest,nest) :: L
double precision, public, DIMENSION(nest,nest) :: U
double precision, public, DIMENSION(nest) :: b
double precision, public :: coeff
double precision, public, DIMENSION(nest) :: d
integer, public :: i
integer, public :: j
integer, public :: k
double precision, public, DIMENSION(nest) :: x

Source Code

 SUBROUTINE inverse(gamma,nest,inv_gamma)

  IMPLICIT NONE 
! Calling variables  
  INTEGER, INTENT(IN) :: nest
  DOUBLE PRECISION, DIMENSION(nest,nest), INTENT(IN) :: gamma 
  DOUBLE PRECISION, DIMENSION(nest,nest), INTENT(OUT) :: inv_gamma
! Local variables  
  DOUBLE PRECISION, DIMENSION(nest,nest)  :: A,L,U
  DOUBLE PRECISION, DIMENSION(nest)  ::  b,d,x
  DOUBLE PRECISION ::  coeff
  INTEGER :: i, j, k
  
! step 0: initialization for matrices 
  A=gamma
  inv_gamma=0.
  L=0.
  U=0.
  b=0.
  d=0.
  x=0.
  coeff=0.
  

! step 1: forward elimination
  DO k=1, nest-1
     DO i=k+1,nest
        coeff=A(i,k)/A(k,k)
        
        L(i,k) = coeff
        DO j=k+1,nest
           A(i,j) = A(i,j)-coeff*A(k,j)
        ENDDO
     ENDDO
  ENDDO
  

! Step 2: prepare L and U matrices 
! L matrix is A matrix of the elimination coefficient
! + the diagonal elements are 1.0
  DO i=1,nest
     L(i,i) = 1.
  ENDDO
! U matrix is the upper triangular part of A
  DO j=1,nest
     DO i=1,j
        U(i,j) = A(i,j)
     ENDDO
  ENDDO

! Step 3: compute columns of the inverse matrix C
  DO k=1,nest
     b(k)=1.0
     d(1) = b(1)
! Step 3a: Solve Ld=b using the forward substitution
     DO i=2,nest
        d(i)=b(i)
        DO j=1,i-1
           d(i) = d(i) - L(i,j)*d(j)
        ENDDO
     ENDDO
! Step 3b: Solve Ux=d using the back substitution
     x(nest)=d(nest)/U(nest,nest)
     DO i = nest-1,1,-1
        x(i) = d(i)
        DO j=nest,i+1,-1
           x(i)=x(i)-U(i,j)*x(j)
        ENDDO
        x(i) = x(i)/u(i,i)
     ENDDO
! Step 3c: fill the solutions x(nest) into column k of C
     DO i=1,nest
        inv_gamma(i,k) = x(i)
     ENDDO
     b(k)=0.0
  ENDDO
 END SUBROUTINE inverse